Based on: J. R. Johansson (robert@riken.jp), http://jrjohansson.github.io, and Eunjong Kim.
In [1]:
from sympy import *
init_printing()
In [2]:
from sympsi import *
from sympsi.boson import *
from sympsi.pauli import *
The Jaynes-Cummings model is one of the most elementary quantum mechanical models light-matter interaction. It describes a single two-level atom that interacts with a single harmonic-oscillator mode of a electromagnetic cavity.
The Hamiltonian for a two-level system in its eigenbasis can be written as
$$ H = \frac{1}{2}\omega_q \sigma_z $$and the Hamiltonian of a quantum harmonic oscillator (cavity) is
$$ H = \hbar\omega_r (a^\dagger a + 1/2) $$The atom interacts with the electromagnetic field produced by the cavity mode $a + a^\dagger$ through its dipole moment. The dipole-transition operators is $\sigma_x$ (which cause a transition from the two dipole states of the atom). The combined atom-cavity Hamiltonian can therefore be written in the form
$$ H = \hbar\omega_r (a^\dagger a + 1/2) + \frac{1}{2}\hbar\Omega\sigma_z + \hbar g\sigma_x(a + a^\dagger) $$Although the Jaynes-Cumming Hamiltonian allow us to evolve the given initial state according to the Schrödinger Equation, in an experiment we would like to predicte the response of the coupled cavity-qubit system under the influence of driving fields for the cavity and qubit, and also account for the effects of dissipation and dephasing (not treated here)
The external coherent-state input may be incorporated in the Jaynes-Cummings Hamiltonian by addition of terms involving the amplitude of the driving field $\vec{E_d} \left(\vec{E_s}\right)$ and it's frequency $\omega_d\left(\omega_s\right)$
$$ H_{cavity} = E_d \left(e^{i\omega_dt}a +e^{-i\omega_dt}\right) $$$$ H_{qubit} = E_s \left(e^{i\omega_st}a +e^{-i\omega_st}\right) $$To obtain the Jaynes-Cumming Hamiltonian
$$ H = \hbar\omega_r (a^\dagger a + 1/2) %-\frac{1}{2}\Delta\sigma_x + \frac{1}{2}\hbar\Omega\sigma_z + \hbar g(\sigma_+ a + \sigma_- a^\dagger) $$we also need to perform a rotating-wave approximation which simplifies the interaction part of the Hamiltonian. In the following we will begin with looking at how these two Hamiltonians are related.
To represent the atom-cavity Hamiltonian in SymPy we creates an instances of the operator classes BosonOp and SigmaX, SigmaY, and SigmaZ, and use these to construct the Hamiltonian (we work in units where $\hbar = 1$).
In [67]:
omega_r, omega_q, g, Delta_d, Delta_s, t, x, chi, Hsym = symbols("omega_r, omega_q, g, Delta_d, Delta_s, t, x, chi, H")
A, B = symbols("A,B") # Electric field amplitude
omega_d, omega_s = symbols("omega_d, omega_s")
Delta = symbols("Delta")
In [4]:
sx, sy, sz, sm, sp = SigmaX(), SigmaY(), SigmaZ(), SigmaMinus(), SigmaPlus()
a = BosonOp("a")
In [59]:
H = omega_r * Dagger(a) * a + omega_q/2 * sz
H_int = g * sx * (a + Dagger(a))
H_drive_r = A * (exp(I*omega_d*t)*a + exp(-I*omega_d*t)*Dagger(a))
H_drive_q = B * (exp(I*omega_s*t)*sm + exp(-I*omega_d*t)*sp)
H_total = H+ H_drive_r + H_drive_q+ H_int
Eq(Hsym, H+ H_drive_r + H_drive_q+ H_int)
Out[59]:
Unitary transformation to interaction picture
In [60]:
U = exp(I * omega_r * t * Dagger(a)*a)
U
Out[60]:
In [61]:
H1 = hamiltonian_transformation(U, H_total.expand())
H1
Out[61]:
In [63]:
U = exp(I * omega_q * t * sp * sm)
U
Out[63]:
In [64]:
H2 = hamiltonian_transformation(U, H1.expand())
H2 = H2.subs(sx, sm + sp).expand()
H2 = powsimp(H2)
H2
Out[64]:
In [65]:
# trick to simplify exponents
def simplify_exp(e):
if isinstance(e, exp):
return exp(simplify(e.exp.expand()))
if isinstance(e, (Add, Mul)):
return type(e)(*(simplify_exp(arg) for arg in e.args))
return e
In [68]:
H3 = simplify_exp(H2).subs(-omega_r + omega_q, Delta)
H3
Out[68]:
Now, in the rotating-wave approximation we can drop the fast oscillating terms containing the factors $e^{\pm i(\omega_q + \omega_r)t}$
In [71]:
H4 = drop_terms_containing(H3, [exp( I * (omega_q + omega_r) * t),
exp(-I * (omega_q + omega_r) * t)])
H4 = drop_c_number_terms(H4.expand())
Eq(Hsym, H4)
Out[71]:
This is the interaction term of in the Jaynes-Cumming model in the interaction picture. If we transform back to the Schrödinger picture we have:
In [72]:
U = exp(-I * omega_r * t * Dagger(a) * a)
H5 = hamiltonian_transformation(U, H4.expand())
H5
Out[72]:
In [73]:
U = exp(-I * omega_q * t * sp * sm)
H6 = hamiltonian_transformation(U, H5.expand())
H6
Out[73]:
In [76]:
H7 = simplify_exp(H6).subs(Delta, omega_q - omega_r)
H7 = simplify_exp(powsimp(H7)).expand()
H7 = drop_c_number_terms(H7)
H = collect(H7, [A,B,g])
Eq(Hsym, H)
Out[76]:
First we apply the unitary transformation $U = e^{i \omega_d a^\dagger a t}$:
In [77]:
U = exp(I * Dagger(a) * a * omega_d * t)
U
Out[77]:
In [78]:
H1 = hamiltonian_transformation(U, H, independent=True)
H1
Out[78]:
We can now perform a rotating-wave approximation (RWA) by eliminating all terms that rotate with frequencies $2\omega_d$:
In [79]:
H2 = drop_terms_containing(H1.expand(), [exp(-2*I*omega_d*t), exp(2*I*omega_d*t)])
Eq(Symbol("H_{rwa}"), H2)
Out[79]:
Introduce the detuning $\Delta_d = \omega_r - \omega_d$:
In [81]:
H3 = H2.subs(omega_r, Delta_d + omega_d).expand()
H3
Out[81]:
Second we apply the unitary transformation $U = e^{i \omega_s \sigma_+ \sigma_- t}$:
In [82]:
U = exp(I * sp*sm* omega_d * t)
U
Out[82]:
In [83]:
H4 = hamiltonian_transformation(U, H3, independent=True)
H4
Out[83]:
In [ ]:
In [ ]:
In [ ]:
Introduce the detuning $\Delta_s = \omega_q - \omega_s$:
In [85]:
H5 = H4.subs(omega_q, Delta_s + omega_s).expand()
H5
Out[85]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
Substitute $ \sigma_x = \sigma_+ + \sigma_-$
In [38]:
H1 = H_total.subs(sx, sm + sp).expand()
H1
Out[38]:
In [39]:
U = exp(I * omega_d * t * (Dagger(a) * a) + I * omega_s * t * (sp * sm))
U.expand()
Out[39]:
In [40]:
H2 = hamiltonian_transformation(U, H1.expand())
H2
Out[40]:
In [41]:
H2a = powsimp(H2)
H2a
Out[41]:
In [42]:
# trick to simplify exponents
def simplify_exp(e):
if isinstance(e, exp):
return exp(simplify(e.exp.expand()))
if isinstance(e, (Add, Mul)):
return type(e)(*(simplify_exp(arg) for arg in e.args))
return e
In [43]:
F = symbols("F")
H2b = simplify_exp(H2a).subs(omega_d+omega_s,F)
H2b
Out[43]:
Now, in the rotating-wave approximation we can drop the fast oscillating terms containing the factors $e^{\pm i(\omega_d + \omega_s)t}$
In [50]:
H2c = H2b.expand()
H2c
H2d= H2c.collect([Dagger(a)*a,sz,g])
H2d
Out[50]:
In [51]:
H2e = H2d.subs(-omega_d+omega_r,Delta_d)
H2e
Out[51]:
In [53]:
H2f = H2e.subs(omega_q/2 - omega_s/2, Delta_s/2)
H2f
Out[53]:
In [55]:
H2g = H2f.expand()
H3 = drop_terms_containing(H2g, [exp( I * (F) * t),
exp(-I * (F) * t)])
H3 = drop_c_number_terms(H3.expand())
Eq(Hsym, H3)
Out[55]:
In [55]:
This is the interaction term of in the Jaynes-Cumming model in the interaction picture. If we transform back to the Schrödinger picture we have:
In [56]:
U = exp(-I * omega_d * t * Dagger(a) * a - I * omega_s * t * sp*sm)
H4 = hamiltonian_transformation(U, H3.expand())
H4
Out[56]:
In [57]:
H4a = H4.expand()
H4a
Out[57]:
In [58]:
H4b = collect(H4a, [Dagger(a)*a,sz,g, A, B])
H4b
Out[58]:
In [30]:
Out[30]:
In [ ]:
In [26]:
H6a = collect(H6, [Dagger(a)*a,sz,g ]);H6a
In [55]:
H6b = powsimp(H6a)
H5 = drop_terms_containing(H6b, [exp( I * (omega_s) * t),
exp(-I * (omega_s) * t)])
H5 = drop_c_number_terms(H5.expand())
H5
Out[55]:
In [57]:
H4 = drop_terms_containing(H5, [t**5, t**4,t**3,t**2,t])
H5 = drop_c_number_terms(H4.expand())
# Eq(Hsym, H5)H
H5
Out[57]:
In [58]:
H = collect(H5, [Dagger(a)*a,sx, g])
Eq(Hsym, H)
Out[58]:
This is the Jaynes-Cumming model give above, and we have now seen that it is obtained to the dipole interaction Hamiltonian through the rotating wave approximation.
In the dispersive regime, where the two-level system is detuned from the cavity by much more than the interaction strength, $\Delta \gg g$, an effective Hamiltonian can be dervied which describes the Stark shift of the two-level system (which depends on the number of photons in the cavity) and the frequency shift of the cavity (which depend on the state of the two-level system).
This effective Hamiltonian, which is correct up to second order in the small paramter $g/\Delta$, is obtained by performing the unitary transformation
$$ U = e^{\frac{g}{\Delta}(a \sigma_- - a^\dagger \sigma_+)} $$
In [20]:
U = exp((x * (a * sp - Dagger(a) * sm)).expand())
U
Out[20]:
In [21]:
#H1 = unitary_transformation(U, H, allinone=True, expansion_search=False, N=3).expand()
#H1 = qsimplify(H1)
#H1
In [22]:
H1 = hamiltonian_transformation(U, H, expansion_search=False, N=3).expand()
H1 = qsimplify(H1)
H1
Out[22]:
In [23]:
H2 = drop_terms_containing(H1.expand(), [x**3, x**4])
H2
Out[23]:
In [24]:
H3 = H2.subs(x, g/Delta)
H3
Out[24]:
In [25]:
H4 = drop_c_number_terms(H3)
H4
Out[25]:
In [26]:
H5 = collect(H4, [Dagger(a) * a, sz])
H5
Out[26]:
Now move to a frame co-rotating with the qubit and oscillator frequencies:
In [27]:
H5.expand()
Out[27]:
In [28]:
U = exp(I * omega_r * t * Dagger(a) * a)
In [29]:
H6 = hamiltonian_transformation(U, H5.expand()); H6
Out[29]:
In [30]:
U = exp(I * Omega * t * Dagger(sm) * sm)
In [31]:
H7 = hamiltonian_transformation(U, H6.expand()); H7
Out[31]:
Now, since we are in the dispersive regime $|\Omega-\omega_r| \gg g$, we can do a rotating-wave approximation and drop all the fast rotating terms in the Hamiltonian above:
In [32]:
H8 = drop_terms_containing(H7, [exp(I * omega_r * t), exp(-I * omega_r * t),
exp(I * Omega * t), exp(-I * Omega * t)])
H8
Out[32]:
In [33]:
H9 = qsimplify(H8)
H9 = collect(H9, [Dagger(a) * a, sz])
H9
Out[33]:
Now move back to the lab frame:
In [34]:
U = exp(-I * omega_r * t * Dagger(a) * a)
In [35]:
H10 = hamiltonian_transformation(U, H9.expand()); H10
Out[35]:
In [36]:
U = exp(-I * Omega * t * Dagger(sm) * sm)
In [37]:
H11 = hamiltonian_transformation(U, H10.expand()); H11
Out[37]:
In [38]:
H12 = qsimplify(H11)
H12 = collect(H12, [Dagger(a) * a, sz])
H12 = H12.subs(omega_r, Omega-Delta).expand().collect([Dagger(a)*a, sz]).subs(Omega-Delta,omega_r)
Eq(Hsym, H12)
Out[38]:
This is the Hamiltonian of the Jaynes-Cummings model in the the dispersive regime. It can be interpreted as the resonator having a qubit-state-dependent frequency shift, or alternatively that the qubit is feeling a resonator-photon-number dependent Stark-shift.
In [35]:
%reload_ext version_information
%version_information sympy, sympsi
Out[35]: